module oct_pet
  !
  ! PET Geometry Definitions
  !
  !
  ! LICENSE AGREEMENT:
  !
  ! You agree that we don't agree to anything.
  !
  !
  ! DISCLAIMER OF WARRANTIES:
  !
  ! This software is provided to you "AS IS," and the authors and
  ! contributors disclaim any warranty or liability obligations to you
  ! of any kind, implied or express, including but not limited to the
  ! warranties of merchantability, fitness for a particular purpose,
  ! and non-infringement.
  ! 
#define green(i) (iand(mana,i)/=0)
  use ground_control
  implicit none
  private
  public setSysmat
  public initOctPet
  public readHist, readList
  public jrt2j
  public getnt
  public jrt2xy
  public x2jx
  public y2jy
  public scaleSysmat
  public writeHist

  !
  real,parameter::pi_4  = atan(1.0_dp)
  real,parameter::pi_2  = pi_4*2.0_dp
  real,parameter::pi    = pi_4*4.0_dp
  real,parameter::pi3_2 = pi_4*6.0_dp
  real,parameter::pi2   = pi_4*8.0_dp

  !
  integer nmaxth

  ! Here we assume crystal width as 2 mm
  real,parameter::          cdt    = 2. 
  real,parameter,public::   cda    = 2.
  real,parameter::          theH   = min(cda,cdt)/2.
  integer,public::          mt     = 0     ! number of crystal in tangential
  integer,public::          ma     = 0     ! number of cyrstal in axial
  integer,parameter,public::md     = 8     ! number of heads

  integer,parameter::   maxdd_def  =  2     ! max difference of detector head
  integer,parameter::   maxcd_def  = -1     ! max difference of crystal
  integer::             maxdd      = maxdd_def
  integer,public::      maxcda
  integer,public::      maxcdt
  integer,public,parameter::npid      = md/2  ! <<< 
  integer,parameter::       ngid_def  =  5
  integer,public::          ngid      =  ngid_def
  real::                       D      = 0.  ! crystal to crystal distance
  real::                    rmax      = 0.  ! maximum radius of FOV 

  !
  integer,public::          nj
  integer,parameter,public::ndt  = 2               ! tangential sampling = 8*ndt
  real,public,parameter::          dx   = cdt/2.        ! width of pixel in x
  real,public,parameter::          dy   = cdt/2.        ! width of pixel in y
  real,parameter::                 dr   = min(dx,dy)/2. ! radial pixel width
  real,public::                    dz   = cda/2.0       ! width of pixel in axial
  integer,public::nr   = 0  ! number of pixel in radial 
  integer,public::nx   = 0  ! number of pixel in x
  integer,public::ny   = 0  ! number of pixel in y
  integer,public::nz   = 0  ! number of pixel in axial

  !
  integer,allocatable::minjs(:),maxjs(:)
  integer,allocatable,public::maxjts(:)

  !
  real,target::mws1(1)=[1.]
  real,target::mtos1(1)=[0.]
  real,target::maos1(1)=[0.]
  real,parameter::sqrt3_5=sqrt(3.0_dp/5.0_dp)
  real,parameter::sqrt_3=sqrt(1.0_dp/3.0_dp)
  real,target::mws4(4)=[1.0_dp/4.0_dp,1.0_dp/4.0_dp,1.0_dp/4.0_dp,1.0_dp/4.0_dp]
  real,target::mtos4(4)=[theH*sqrt_3,-theH*sqrt_3,-theH*sqrt_3,theH*sqrt_3]
  real,target::maos4(4)=[theH*sqrt_3, theH*sqrt_3,-theH*sqrt_3,-theH*sqrt_3]

  real,target::mws9(9)=[25.0_dp/324.0_dp,25.0_dp/324.0_dp,25.0_dp/324.0_dp,25.0_dp/324.0_dp,&
       10.0_dp/81.0_dp,10.0_dp/81.0_dp,10.0_dp/81.0_dp,10.0_dp/81.0_dp,&
       16.0_dp/81.0_dp]
  real,target::mtos9(9)=[theH*sqrt3_5,-theH*sqrt3_5,-theH*sqrt3_5,theH*sqrt3_5,&
       theH*sqrt3_5,0.,-theH*sqrt3_5,0.,&
       0.]
  real,target::maos9(9)=[theH*sqrt3_5,theH*sqrt3_5,-theH*sqrt3_5,-theH*sqrt3_5,&
       0.,theH*sqrt3_5,0.,-theH*sqrt3_5,&
       0.]
  real,pointer::mtos(:)
  real,pointer::maos(:)
  real,pointer::mws(:)
  
  public t_sysmat
  public t_block_pair
  public t_grouped_pair

  type t_sysmat
     sequence
     integer n               ! number of voxels
     integer p               ! pointer to jr,jt,jz,jw
  end type t_sysmat
  type(t_sysmat),target,allocatable,public::sysmat(:,:,:,:)
  integer ntsme,nmaxsme,szsysmat ! number of total/max sysmat elem, size of sysmat
  real maxsmw, minsmw
  integer,parameter::size_of_sysmat_element=4+4+8 ! j, jz, w
  
  type t_block_pair
     sequence
     integer id1
     integer id2
     integer rot
  end type t_block_pair
  type t_grouped_pair
     sequence
     integer pid
     integer gid
  end type t_grouped_pair
  type(t_block_pair),public,allocatable::idt(:,:)
  type(t_grouped_pair),public,allocatable::ridt(:,:)
  integer,public::numidt

  real,allocatable,target::tmp_stat(:,:)
  integer nbufssm

  integer lismf ! line integral sampling multiplication factor
  integer nmr   ! number of multi ray

contains

  subroutine echo(mess,no_nl)
    use suty
    character*(*),intent(in)::mess
    integer,intent(in),optional::no_nl
    character(MAXPROCLEN),save::theProc=""
    character(MAXPROCLEN) proc
    if(.not.green(OPT_VERBOSE)) return
    proc=getProcname(mess)
    if(proc=="") then
       proc=theProc
    else
       theProc=proc
       proc=""
    end if
    call echo_green(mess,proc,"op",.false.,no_nl)
  end subroutine echo

  subroutine initOctPet(penv)
    use omp_lib
    use suty
    type(t_opet_env),intent(in)::penv

    nmaxth=omp_get_max_threads()
    call echo("initOctPet: number of threads = "//trim(i2a(nmaxth)))

    ma=penv%na
    mt=penv%nt
    call echo("Crystal array: "//trim(i2a(mt))//" x "//trim(i2a(ma))&
         //" ("//trim(f2a(cdt,-4))//" x "//trim(f2a(cda,-4))//" mm)")
    D=penv%D
    rmax = D/2.0_dp - 2.0_dp

    call echo("Detector separation = "//trim(f2a(penv%D,-4)))
    if(penv%D<cdt*mt*(1.+sqrt(2.))) call echo_red("initOctPet: FOV too small")
    
    maxdd = maxdd_def
    if(penv%maxdd>0) maxdd=penv%maxdd
    if(penv%maxcd>0) then
       maxcda=min(penv%maxcd,ma-1)
       maxcdt=min(penv%maxcd,mt-1)       
    else
       maxcda=ma-1
       maxcdt=mt-1
    end if

    if(penv%R/=0) then
       if(penv%R>rmax) call echo_red("**** FOV radius too large: "//trim(f2a(penv%R,-4)))
       if(penv%R<=dr) call echo_red("*** FOV radius too small: "//trim(f2a(penv%R,-4)))
       rmax = penv%R
    end if

    if(penv%dz>0) dz = penv%dz
    nr   = floor(rmax/dr)      
    rmax = dr*nr
    call echo("FOV radius = "//trim(f2a(rmax,-4)))
    nx   = nint(real(rmax,dp)/dx)*2+1
    ny   = nint(real(rmax,dp)/dy)*2+1 
    nz   = max(1,ma*nint(cda/dz))+1

    if(maxdd>2) maxdd=2
    ngid = maxdd*2+1
    call echo("maxcda, maxcdt, maxdd = "//trim(i2a(maxcda))//", "//trim(i2a(maxcdt))//", "//trim(i2a(maxdd)))
        
    call setIdt

    call setRidt

    call setProjGeom
    call echo("cylindrical sampling  = "//trim(i2a(nj)))
    call echo("radial sanpling = "//trim(i2a(nr))//" ("//trim(f2a(dr,-4))//" mm)")
    call echo("image dimension: "//trim(i2a(nx))//" x "//trim(i2a(ny))//" x "//trim(i2a(nz)))
    call echo("pixel size: "//trim(f2a(dx,-4))//" x "//trim(f2a(dy,-4))//" mm")
    call echo("slice thickness = "//trim(f2a(dz,-4))//" mm")

  contains
    
    subroutine setIdt

      if(allocated(idt)) deallocate(idt)
      allocate(idt(npid,ngid))
      numidt=npid*ngid

      idt%id1=0
      idt%id2=0
      idt%rot=0
      
      ! group 1, dd=0
      idt(1,1)%id1=1; idt(2,1)%id1=2; idt(3,1)%id1=3; idt(4,1)%id1=4
      idt(1,1)%id2=5; idt(2,1)%id2=6; idt(3,1)%id2=7; idt(4,1)%id2=8
      idt(1,1)%rot=0; idt(2,1)%rot=1; idt(3,1)%rot=2; idt(4,1)%rot=3

      if(maxdd==0) return
      
      ! group 2, dd=1
      idt(1,2)%id1=1; idt(2,2)%id1=3; idt(3,2)%id1=5; idt(4,2)%id1=-2
      idt(1,2)%id2=4; idt(2,2)%id2=6; idt(3,2)%id2=8; idt(4,2)%id2=7
      idt(1,2)%rot=0; idt(2,2)%rot=2; idt(3,2)%rot=4; idt(4,2)%rot=6
      
      ! group 3 dd=-1
      idt(1,3)%id1=2; idt(2,3)%id1=4; idt(3,3)%id1=-1; idt(4,3)%id1=-3
      idt(1,3)%id2=5; idt(2,3)%id2=7; idt(3,3)%id2=6 ; idt(4,3)%id2=8  
      idt(1,3)%rot=1; idt(2,3)%rot=3; idt(3,3)%rot=5 ; idt(4,3)%rot=7  
      
      if(maxdd==1) return

      ! group 4, dd=2
      idt(1,4)%id1=1; idt(2,4)%id1=3; idt(3,4)%id1=5; idt(4,4)%id1=-1
      idt(1,4)%id2=3; idt(2,4)%id2=5; idt(3,4)%id2=7; idt(4,4)%id2=7 
      idt(1,4)%rot=0; idt(2,4)%rot=2; idt(3,4)%rot=4; idt(4,4)%rot=6 
            
      ! group 5, dd=-2
      idt(1,5)%id1=2; idt(2,5)%id1=4; idt(3,5)%id1=6; idt(4,5)%id1=-2
      idt(1,5)%id2=4; idt(2,5)%id2=6; idt(3,5)%id2=8; idt(4,5)%id2=8 
      idt(1,5)%rot=1; idt(2,5)%rot=3; idt(3,5)%rot=5; idt(4,5)%rot=7 
      
    end subroutine setIdt
    
    subroutine setRidt
      integer ip,ig
      integer id1,iip
      type(t_block_pair) bp
      pointer(pbp,bp)
      if(allocated(ridt)) deallocate(ridt)
      allocate(ridt(md,md))
      ridt%pid=0
      ridt%gid=0
      do ig=1,ngid
         do ip=1,npid
            pbp=loc(idt(ip,ig))
            id1=idt(ip,ig)%id1
            iip=ip
            if(id1<0) then
               id1=-id1
               iip=-iip
            end if 
            ridt(id1,bp%id2)%pid=iip
            ridt(id1,bp%id2)%gid=ig           
         end do
      end do
    end subroutine setRidt

    subroutine setProjGeom
      integer j,jr
      
      nj=jrt2j(nr,getnt(nr))
      if(allocated(minjs)) deallocate(minjs)
      if(allocated(maxjs)) deallocate(maxjs)
      if(allocated(maxjts)) deallocate(maxjts)
      allocate(minjs(0:nr),maxjs(0:nr),maxjts(0:nr))
      
      minjs(0)=0
      maxjs(0)=1
      maxjts(0)=1
      j=1
      do jr=1,nr
         minjs(jr)=j
         maxjts(jr)=getnt(jr)
         j=j+maxjts(jr)
         maxjs(jr)=j-1
      end do
      
    end subroutine setProjGeom

  end subroutine initOctPet

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!! Sysmat I/O funtions
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  integer function loadSysmat(f)
    character*(*),intent(in)::f
    integer istat
    loadSysmat=-2
    open(unit=1,file=f,status="old",access="stream",iostat=istat)
    if(istat/=0) return
    loadSysmat=-1
    if(checkSysmatHdr()) then
       call initSysmat
       call readSysmatElement
    end if
    close(1)
    if(istat/=0) return
    loadSysmat=0
    contains
      logical function checkSysmatHdr()
        integer*4 nx_,ny_,nr_,nz_,ndt_
        real*4 dr_        
        checkSysmatHdr=.false.
        read(1,iostat=istat) nx_,ny_,nr_,nz_,ndt_,dr_
        if(istat/=0) return
        if(nx_==nx.and.ny_==ny.and.&
             nr_==nr.and.nz_==nz.and.ndt_==ndt.and.dr_==real(dr,kind=4)) checkSysmatHdr=.true.
      end function checkSysmatHdr
      subroutine readSysmatElement
        integer*1 b(*)
        integer s
        integer id,it1,it2,ia2
        type(t_sysmat),pointer::sm
        pointer(p,b)
        do id=1,maxdd+1
           do ia2=1,ma
              do it2=1,mt
                 do it1=1,mt
                    sm => sysmat(it1,it2,ia2,id)
                    read(1,iostat=istat) sm%n
                    sm%p=0
                    if(istat/=0) return
                    if(sm%n/=0) then
                       ntsme=ntsme+sm%n
                       nmaxsme=max(nmaxsme,sm%n)
                       s=sm%n*size_of_sysmat_element
                       szsysmat=szsysmat+s
                       p=malloc(s)
                       sm%p=p
                       read(1,iostat=istat) b(1:s)
                       if(istat/=0) exit
                    end if
                 end do
              end do
           end do
        end do
      end subroutine readSysmatElement      
  end function loadSysmat

  integer function saveSysmat(f)
    character*(*),intent(in)::f
    integer istat
    saveSysmat=-2
    open(unit=1,file=f,status="replace",access="stream",iostat=istat)
    if(istat/=0) return
    saveSysmat=-1
    call writeSysmatHdr
    if(istat==0) call writeSysmatElement
    close(1)
    if(istat/=0) return
    saveSysmat=0
  contains
    subroutine writeSysmatHdr
      write(1,iostat=istat) int(nx,kind=4),int(ny,kind=4),&
           int(nr,kind=4),int(nz,kind=4),int(ndt,kind=4),real(dr,kind=4)
    end subroutine writeSysmatHdr
    subroutine writeSysmatElement
      integer id,it1,it2,ia2
      integer*1 b(*)
      type(t_sysmat),pointer::sm
      pointer(p,b)    
      do id=1,maxdd+1
         do ia2=1,ma
            do it2=1,mt
               do it1=1,mt
                  sm => sysmat(it1,it2,ia2,id)
                  write(1,iostat=istat) int(sm%n,kind=4)
                  if(istat/=0) exit
                  if(sm%p/=0) then
                     p=sm%p
                     write(1,iostat=istat) b(1:sm%n*size_of_sysmat_element)
                     if(istat/=0) exit
                  end if
               end do
            end do
         end do
      end do
    end subroutine writeSysmatElement
  end function saveSysmat
  
  subroutine deallocSysmat()
    integer id,it1,it2,ia2
    type(t_sysmat),pointer::sm
    if(.not.allocated(sysmat)) return
    do id=1,maxdd+1
       do ia2=1,ma
          do it2=1,mt
             do it1=1,mt
                sm => sysmat(it1,it2,ia2,id)
                if(sm%p==0) cycle
                call free(sm%p)
             end do
          end do
       end do
    end do
    deallocate(sysmat)
  end subroutine deallocSysmat

  subroutine initSysmat()
    integer id,it1,it2,ia2
    type(t_sysmat),pointer::sm
    call deallocSysmat
    allocate(sysmat(nx,nx,ny,maxdd+1))
    do id=1,maxdd+1
       do ia2=1,ma
          do it2=1,mt
             do it1=1,mt
                sm => sysmat(it1,it2,ia2,id)
                sm%p=0
                sm%n=0
             end do
          end do
       end do
    end do
    ntsme=0
    nmaxsme=0
    szsysmat=0
  end subroutine initSysmat
  
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!! Sysmat evaluator
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  
  subroutine getCrystalPos(id,it,ia,im,xs,vo)
    use suty
    integer,intent(in)::id,it,ia,im
    real,intent(out),target::xs(3)
    real,intent(in),optional::vo
    real,pointer::x,y,z
    real::xx=0,yy=0
    x => xs(1); y => xs(2); z => xs(3)
    !
    x=it*cdt-mt*cdt/2.0_dp-cdt/2.0_dp+mtos(im)
    y=D/2
    if(present(vo)) y=y+vo
    z=ia*cda-ma*cda/2.0_dp-cda/2.0_dp+maos(im)
    !
    if(id==1) return
    select case(id)
    case(1)
       return
    case(2)
       xx=( x - y)/sqrt(2.0_dp)
       yy=(+x + y)/sqrt(2.0_dp)
    case(3)
       xx = -y
       yy = +x
    case(4)
       xx=(-x - y)/sqrt(2.0_dp)
       yy=(+x - y)/sqrt(2.0_dp)
    case(5)
       xx=-x
       yy=-y
    case default
       call echo_red("getCrystalPos: invalid id = "//trim(i2a(id)))
    end select
    x=xx
    y=yy
  end subroutine getCrystalPos
  
  !
  subroutine scaleSysmat(a)
    use suty
    real,intent(in)::a
    integer it1,it2,ia2,idd
    type(t_sysmat),pointer:: sm
    real w(*)
    pointer(p,w)
    call echo("scaleSysmat: scaling factor = "//trim(f2a(a)))
    do idd=1,maxdd+1
       do ia2=1,ma
          do it2=1,mt
             do it1=1,mt
                sm => sysmat(it1,it2,ia2,idd)             
                if(sm%n==0) cycle
                p = sm%p+sm%n*(4+4)
                w(1:sm%n)=w(1:sm%n)*a
             end do
          end do
       end do
    end do
  end subroutine scaleSysmat

  !
  subroutine setSysmat(n,m)
    use omp_lib
    use suty
    integer,intent(in)::n,m    
    integer ia2,ia2_end
    nmr = n
    lismf = m
    call echo("setSysmat: number of multi-ray = "//trim(i2a(nmr)))
    call echo("number of line integral sampling = "//trim(i2a(lismf)))
    select case(nmr)
    case(4)
       mtos => mtos4
       maos => maos4
       mws  =>  mws4
    case(9)
       mtos => mtos9
       maos => maos9
       mws  =>  mws9
    case(1)
       mtos => mtos1
       maos => maos1
       mws  =>  mws1      
    case default
       call echo_red("setSysmat: invalid nmr.")
    end select
    call echo("Begin.",1)
    call initSysmat
    if(allocated(tmp_stat)) deallocate(tmp_stat)
    allocate(tmp_stat(6,nmaxth))
    tmp_stat=0
    ia2_end=ma-(ma-maxcda)+1
    call doSetSysmat(ia2_end,0,nmr,fake=.true.)
    nbufssm = int(tmp_stat(6,1))
    tmp_stat=0
    tmp_stat(5,:)=huge(0.0_dp)
    tmp_stat(4,:)=-huge(0.0_dp)
    !$omp parallel do schedule(dynamic) private(ia2)
    do ia2=1,ia2_end
       call doSetSysmat(ia2,0,nmr)
       call doSetSysmat(ia2,1,nmr)
       call doSetSysmat(ia2,2,nmr)
    end do
    !$omp end parallel do
    call echo("done.",2)
    ntsme=int(sum(tmp_stat(1,:)))
    nmaxsme=int(sum(tmp_stat(2,:)))
    szsysmat=int(sum(tmp_stat(3,:)))
    minsmw=minval(tmp_stat(5,:))
    maxsmw=maxval(tmp_stat(4,:))    
    call echo("total number of element:   "//trim(i2a(ntsme,12)))
    call echo("maximum number of element: "//trim(i2a(nmaxsme,12)))
    call echo("size of system matrix:     "//trim(i2a(szsysmat,12))//" Byte")
    call echo("maximum buffer number used:"//trim(i2a(int(maxval(tmp_stat(6,:))),12)))
    call echo("maximum value of sysmat:   "//trim(f2a(maxsmw)))
    call echo("minimum value of sysmat:   "//trim(f2a(minsmw)))
  end subroutine setSysmat
  !
  subroutine doSetSysmat(ia2,dd,nmr,fake)
    use omp_lib
    use suty
    integer,intent(in)::ia2
    integer,intent(in)::dd
    integer,intent(in)::nmr
    logical,intent(in),optional::fake
    real,pointer::stat(:)
    integer ith
    integer idd
    integer id1,id2
    integer it1,ia1,it2
    real x1s(3),x2s(3)
    real xr1s(3),xr2s(3)
    real xs(3),xvs(3)
    real M(3),L(3)
    real nL,theL,norm
    real k,dk
    integer nk    
    integer j,jz,jr,jt
    real t,r
    real w
    integer im1,im2
    integer ie,n,p
    integer ntot
    integer nmax
    integer nbufmax
    integer rsz
    real minw,maxw
    real d2
    type t_jbufel
       sequence
       integer*4 j
       integer*4 jz
    end type t_jbufel
    type(t_sysmat),pointer::sm
    type(t_jbufel),allocatable::jbuf(:)
    real,allocatable:: wbuf(:)
    integer*8 jjz(*)
    pointer(pjjz,jjz)
    logical run
    if(dd>maxdd) return
    ith=omp_get_thread_num()+1
    stat => tmp_stat(:,ith)
    ntot=0
    nmax=0
    rsz =0
    nbufmax=0
    minw = huge(1.0_dp)
    maxw = -1.0_dp
    run=.not.present(fake).or..not.fake
    if(run) allocate(jbuf(nbufssm),wbuf(nbufssm))   
    pjjz=loc(jbuf)
    idd=dd+1
    if(run.and.green(OPT_VERBOSE)) then
       write(*,10) "."
10     format(a,$)
    end if
    !------------------------------------
    id1=1
    id2=5-dd
    !------------------------------------
    ia1=1
    tangential_1: do it2=1,mt
       tangential_2: do it1=1,mt
          if(abs(it1-it2)>maxcdt) cycle
          ie=0
          sm => sysmat(it1,it2,ia2,idd)
          multiray_2: do im2=1,nmr
             call getCrystalPos(id2,it2,ia2,im2,xr2s)
             x2s=xr2s
             multiray_1: do im1=1,nmr
                call getCrystalPos(id1,it1,ia1,im1,xr1s)
                x1s=xr1s
                M=(x1s+x2s)/2.0_dp
                L=(x1s-x2s)/2.0_dp
                theL=sqrt(sum(L**2.0_dp))
                nL=theL*lismf
                nk=2*nint(nL/dr)+1
                dk=2.0_dp/real(nk,dp)
                k=-1
                line_integral: do
                   k=k+dk
                   if(k>=1) exit
                   xs=L*k+M
                   call xy2rt(xs(1),xs(2),r,t) ! the nearest neighbor
                   if(r>rmax) cycle
                   call z2jz(xs(3),jz) ! same as above
                   if(jz<=0.or.jz>nz) then
                      call warn("jz under/over-flow",jz)
                      cycle
                   end if
                   ! need some kind of weight here (theL/2.d0*(1+abs(k)))**(-2)!*sqrt(r+dr*4)
                   norm=1 
                   xvs(3)=jz2z(jz)
                   jr=r2jr(r)
                   jt=t2jt(t,jr)
                   call jrt2xy(jr,jt,xvs(1),xvs(2))
                   d2=sum((xvs-xs)**2)
                   call saveOnBuf
                   !
                   if(green(OPT_USE_SNN)) then
                      if(d2/=0.and.findNextNearest(jr,jt,jz,xs,d2)) call saveOnBuf
                   end if
                end do line_integral
             end do multiray_1
          end do multiray_2
          
          if(.not.run) cycle
          n=reduced(ie)
          call store
          ntot=ntot+n
          nmax=max(n,nmax)
          rsz=rsz+n*size_of_sysmat_element
          
       end do tangential_2
    end do tangential_1

    stat(1)=ntot+stat(1)
    stat(2)=nmax+stat(2)
    stat(3)=rsz+stat(3)
    stat(4)=max(maxw,stat(4))
    stat(5)=min(minw,stat(5))
    stat(6)=nbufmax
    if(run) deallocate(jbuf,wbuf)

  contains   

    integer function reduced(num_element)
      integer,intent(in)::num_element
      integer ii
      integer kk
      integer jold
      integer per(num_element)
      real wwbuf(num_element)
      call quicksort(jjz(1:num_element),per(1:num_element))
      wwbuf(1:num_element)=wbuf(1:num_element)
      do ii=1,num_element
         wbuf(ii)=wwbuf(per(ii))
      end do
      kk=0
      jold=0
      do ii=1,num_element
         if(jjz(ii)/=jold) then
            kk=kk+1
            if(kk/=ii) then
               jjz(kk)=jjz(ii)
               wbuf(kk)=wbuf(ii)
            end if
            jold=jjz(ii)
         else
            wbuf(kk)=wbuf(kk)+wbuf(ii)
         end if
      end do
      reduced=kk
      maxw=max(maxw,maxval(wbuf(1:kk)))
      minw=min(minw,minval(wbuf(1:kk)))      
    end function reduced
    subroutine swap(a,b)
      real,intent(inout)::a,b
      real c
      c=a
      a=b
      b=c
    end subroutine swap
    subroutine saveOnBuf
      j=jrt2j(jr,jt)
      if(j<=0.or.j>nj) then
         call warn("j under/over-flow",j)
         return
      end if
      ie=ie+1
      if(ie>nbufmax) nbufmax=ie
      if(.not.run) return
      if(ie>nbufssm) &
           call echo_red("Internal buffer overflow nbuf="//trim(i2a(ie)))
!      w=exp(-d2/smw**2.d0)*mws(im1)*mws(im2)
      w=exp(-d2)*mws(im1)*mws(im2)
      wbuf(ie)=w*norm  ! need some kind of weight here
      jbuf(ie)%jz=int(jz,kind=4)
      jbuf(ie)%j=int(j,kind=4)
    end subroutine saveOnBuf
    subroutine store      
      use suty
      integer ii
      integer*4 js(n),jzs(n)
      sm%n=n
      if(n==0) return
      do ii=1,n
         js(ii)=jbuf(ii)%j
         jzs(ii)=jbuf(ii)%jz
      end do
      p=malloc(n*size_of_sysmat_element)
      sm%p=p
      call cpmems4(p,     loc(js),   n)
      call cpmems4(p+n*4, loc(jzs),  n)
      call cpmems8(p+n*8, loc(wbuf), n)            
    end subroutine store
    subroutine warn(mess,indx)
      character*(*),intent(in)::mess
      integer,intent(in)::indx
      if(.not.green(OPT_DEBUG)) return
      write(*,*) "doSetSysmat: "//trim(mess)//": ",indx
      write(*,*) "at   r1 = ",x1s
      write(*,*) "and  r2 = ",x2s
      write(*,*) "with rk = ",xs
    end subroutine warn
    logical function findNextNearest(jr0,jt0,jz0,x0,dist2)
      integer,intent(inout)::jr0,jt0,jz0
      real,intent(in)::x0(3)
      real,intent(out)::dist2
      real xn(3)
      integer jrs(6)
      integer jts(6)
      integer jzs(6)
      integer ip,ipmin,ip1
      real distmin
      jrs=jr0
      jts=jt0
      jzs=jz0
      ! radial
      if(jr0>0) then
         jrs(1)=jrs(1)+1
         jts(1)=nint(real(jt0,dp)*(1.0_dp+1.0_dp/real(jr0,dp)))
         jrs(2)=jrs(2)-1
         jts(2)=nint(real(jt0,dp)*(1.0_dp-1.0_dp/real(jr0,dp)))
         ! tangential
         if(jt0<8*ndt*jr0) then
            jts(3)=jts(3)+1
         else
            jts(3)=1
         end if
         if(jt0>1) then
            jts(4)=jts(4)-1
         else
            jts(4)=8*ndt*jr0
         end if
         ip1=5
      else
         ip1=1
      end if
      ! axial
      jzs(5)=jzs(5)+1
      jzs(6)=jzs(6)-1
      ipmin=0
      distmin=huge(0.0_dp)
      do ip=ip1,6
         if(jrs(ip)>nr.or.jrs(ip)<=0&
              .or.jzs(ip)>nz.or.jzs(ip)<=0) cycle
         call jrt2xy(jrs(ip),jts(ip),xn(1),xn(2))
         xn(3)=jz2z(jzs(ip))
         dist2=sum((xn-x0)**2)
         if(dist2<distmin) then
            distmin=dist2
            ipmin=ip
         end if
      end do
      if(ipmin/=0) then
         jr0=jrs(ipmin)
         jt0=jts(ipmin)
         jz0=jzs(ipmin)
         dist2=distmin
         findNextNearest=.true.
      else
         findNextNearest=.false.
      end if
    end function findNextNearest
  end subroutine doSetSysmat

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!! Coordinate Transformations and others
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    
  real function getdt(jr)
    integer,intent(in)::jr
    if(jr==0) then
       getdt=0
    else
       getdt=pi_4/(jr*ndt)
    end if
  end function getdt

  integer function getnt(jr)
    integer,intent(in)::jr
    if(jr==0) then
       getnt=1
    else
       getnt=jr*8*ndt
    end if
  end function getnt

  real function getdl(jr)
    integer,intent(in)::jr
    getdl=getdt(jr)*jr2r(jr)
  end function getdl

  subroutine xy2rt(x,y,r,t)
    real,intent(in)::x,y
    real,intent(out)::r,t
    r=sqrt(x**2.0_dp+y**2.0_dp)
    t=atan(y/x)
    if(x<0) then
       t=t+pi
    else if(y<0) then
       t=pi2+t
    end if
  end subroutine xy2rt

  subroutine rt2xy(r,t,x,y)
    real,intent(in)::r,t
    real,intent(out)::x,y
    x=r*cos(t)
    y=r*sin(t)
  end subroutine rt2xy
  subroutine jrt2xy(jr,jt,x,y)
    integer,intent(in)::jr,jt
    real,intent(out)::x,y
    real r,t
    if(jr==0) then
       x=0
       y=0
    else
       r=jr2r(jr)
       t=jrt2t(jr,jt)
       x=r*cos(t)
       y=r*sin(t)
    end if
  end subroutine jrt2xy

  real function jr2r(jr)
    integer,intent(in)::jr
    jr2r=dr*jr
  end function jr2r

  integer function r2jr(r,fr)
    real,intent(in)::r
    real,intent(out),optional::fr
    real v
    v=r/dr
    if(.not.present(fr)) then
       r2jr=nint(v)
    else
       r2jr=floor(v)
       fr=v-r2jr
    end if
  end function r2jr

  subroutine j2jrt(j,jr,jt)
    integer,intent(in)::j
    integer,intent(out)::jr,jt
    if(j==1) then 
       jr=0
       jt=1
    else
       jr=j2jr(j)
       jt=j2jt(j,jr)
    end if
  end subroutine j2jrt

  integer function j2jt(j,jr)
    integer,intent(in)::j,jr
    if(jr>nr) call echo_red( "j2jt: jr>nr")
    j2jt=j-minjs(jr)
  end function j2jt
    
  integer function j2jr(j)
    integer,intent(in)::j
    integer jr
    do jr=1,nr
       if(minjs(jr)>=j) then
          j2jr=jr-1
          return
       end if
    end do
    j2jr=nr
  end function j2jr

  integer function t2jt(t,jr,ft)
    real,intent(in)::t
    integer,intent(in)::jr
    real,intent(out),optional::ft
    real v,dt
    if(jr==0) then
       t2jt=1
       if(present(ft)) ft=0
    else
       dt=pi_4/real(ndt*jr,dp)
       v=t/dt
       if(present(ft)) then
          t2jt=floor(v)
          ft=v-t2jt
       else
          t2jt=nint(v)
       end if
       t2jt=t2jt+1
       if(t2jt>jr*8*ndt) t2jt=1  
    end if
  end function t2jt

  real function jrt2t(jr,jt)
    integer,intent(in)::jr,jt
    ! jt=1:8*jr
    if(jr==0.or.jt==1) then 
       jrt2t=0
    else
       jrt2t=pi_4*real(jt-1,dp)/real(jr*ndt,dp)
    end if
  end function jrt2t
  
  integer function xy2j(x,y)
    real,intent(in)::x,y
    real r,t
    call xy2rt(x,y,r,t)
    xy2j=rt2j(r,t)
  end function xy2j

  integer function rt2j(r,t)
    real,intent(in)::r,t
    integer jr
    jr=r2jr(r)
    rt2j=jrt2j(jr,t2jt(t,jr))
  end function rt2j

  integer function jrt2j(jr,jt)
    integer,intent(in)::jr,jt
    ! jr=0:nr
    ! jt=1:8*ndt*jr
    if(jr==0) then
       jrt2j=1 ! r=t=0
    else if(jt>8*jr*ndt.or.jt<0.or.jr<0) then
       write(*,*) "jr=",jr," jt=",jt," jtmax=",8*jr*ndt
       call echo_red("jrt2j: invalid jr,jt")
       jrt2j=-1
    else
       jrt2j=4*ndt*jr*(jr-1)+jt+1
    end if
  end function jrt2j

  real function jx2x(jx)
    integer,intent(in)::jx
    jx2x=dx*(jx-(nx+1)/2)
  end function jx2x
  real function jy2y(jy)
    integer,intent(in)::jy
    jy2y=dy*(jy-(ny+1)/2)
  end function jy2y
  real function jz2z(jz)
    integer,intent(in)::jz
    jz2z=dz*(jz-(nz+1)/2)
  end function jz2z

  subroutine x2jx(x,jx,fx)
    real,intent(in)::x
    integer,intent(out)::jx
    real,intent(out),optional::fx
    real v
    v=x/dx+(nx+1)/2
    if(.not.present(fx)) then
       jx=nint(v)
    else
       jx=floor(v)
       fx=v-jx
    end if
  end subroutine x2jx
  subroutine y2jy(y,jy,fy)
    real,intent(in)::y
    integer,intent(out)::jy
    real,intent(out),optional::fy
    real v
    v=y/dy+(ny+1)/2
    if(.not.present(fy)) then
       jy=nint(v)
    else
       jy=floor(v)
       fy=v-jy
    end if
  end subroutine y2jy
  subroutine z2jz(z,jz,fz)
    real,intent(in)::z
    integer,intent(out)::jz
    real,intent(out),optional::fz
    real v
    v=z/dz+real(nz+1,dp)/2.0_dp
    if(.not.present(fz)) then
       jz=nint(v)
    else
       jz=floor(v)
       fz=v-jz
    end if
  end subroutine z2jz

  !
  subroutine readHist(f,h_i)
    use suty
    character*(*),intent(in)::f
    real,intent(out)::h_i(:,:,:,:,:,:)
    integer*2 buf(mt,ma,mt,ma,4,ngid)
    integer*4 buf4(mt,ma,mt,ma,4,ngid)
    integer istat

    call echo("readHist: f = "//trim(f))
    open(unit=1,file=f,status="old",access="stream",iostat=istat)
    if(istat/=0) call echo_red("*** cannot open histogram")
    
    if(green(OPT_READ_HIST_INT4)) then      
       read(1,iostat=istat) buf4
       buf=int(buf4,2)
    else
       read(1,iostat=istat) buf
    end if
    if(istat/=0) call echo_red("*** cannot read histogram")
    close(1)
   
    call echo("max/sum = "//trim(i2a(maxval(buf)))//" / "//trim(i2a(sum(buf))))

    h_i=h_i+buf
  end subroutine readHist

  subroutine readList(f,h_i,r_i)
    use suty
    character*(*),intent(in)::f
    real,intent(inout)::h_i(:,:,:,:,:,:)
    real,intent(inout),optional::r_i(:,:,:,:,:,:)
    integer istat
    integer*4 dws(2)
    logical prompt
    integer t1,a1,d1
    integer t2,a2,d2
    integer nread,nacc,nvalid
    integer ip,ig

    call echo("readList: f = "//trim(f))
    open(unit=1,file=f,status="old",access="stream",iostat=istat)
    if(istat/=0) call echo_red("*** cannot open LMF data")

    nread=0
    nvalid=0
    nacc=0
    read: do while(istat==0)
       read(1,iostat=istat) dws(1)
       if(istat/=0) exit
       nread = nread + 4

       if(btest(dws(1),31)) cycle
       if(.not.btest(dws(1),30)) cycle

       isup: do 
          read(1,iostat=istat) dws(2)
          if(istat/=0) exit
          nread = nread + 4

          if(btest(dws(2),31)) cycle read
          if(btest(dws(2),30)) then
             dws(1)=dws(2)
             cycle isup
          else 
             exit isup
          end if
          
          nvalid = nvalid + 1

          prompt=btest(dws(2),28)

          t1=ibits(dws(1),0,6) +1
          a1=ibits(dws(1),6,6) +1
          d1=ibits(dws(1),12,4)+1
          
          if(t1>mt.or.a1>ma.or.d1>md) cycle

          t2=ibits(dws(2),0,6) +1
          a2=ibits(dws(2),6,6) +1
          d2=ibits(dws(2),12,4)+1

          if(t2>mt.or.a2>ma.or.d2>md) cycle

#ifdef _DEBUG_
          if(green(OPT_LMF_FLIP_PMTID)) then
             if(d1/=1) d1=md-d1+2
             if(d2/=1) d2=md-d2+2
          end if

          if(green(OPT_LMF_SWAP_AT)) then
             call swap(a1,t1)
             call swap(a2,t2)
          end if

          if(green(OPT_LMF_FLIP_AXIAL)) then
             a1=ma-a1+1
             a2=ma-a2+1
          end if
          
          if(green(OPT_LMF_FLIP_TAXIAL)) then
             t1=mt-t1+1
             t2=mt-t2+1
          end if
#endif

          if(d1>d2) then
             call swap(a1,a2)
             call swap(t1,t2)
             call swap(d1,d2)
          end if

          ip=ridt(d1,d2)%pid
          ig=ridt(d1,d2)%gid
          if(ip==0.or.ig==0) cycle
          if(ip<0) then
             call swap(a1,a2)
             call swap(t1,t2)
             ip=-ip
          end if
          nacc = nacc +1

          if(prompt) then
             h_i(t1,a1,t2,a2,ip,ig)=h_i(t1,a1,t2,a2,ip,ig)+1
          else if(present(r_i)) then
             r_i(t1,a1,t2,a2,ip,ig)=r_i(t1,a1,t2,a2,ip,ig)+1            
          end if

       end do isup
    end do read
    close(1)

    call echo("bytes read: "//i2a(nread))
    call echo("valid/accepted = "//i2a(nvalid)&
         //" / "//i2a(nacc))
    call echo("Prompt max/sum = "//trim(i2a(int(maxval(h_i))))//" / "//trim(i2a(int(sum(h_i)))))
    if(present(r_i)) call echo("Delayed max/sum = "//trim(i2a(int(maxval(r_i))))//" / "//trim(i2a(int(sum(r_i)))))
  end subroutine readList

  subroutine writeHist(h_i,f)
    real,intent(in)::h_i(:,:,:,:,:,:)
    character*(*),intent(in)::f
    integer*2,allocatable::buf(:,:,:,:,:,:)
    character(MAXPATHLEN) fout
    integer istat
    allocate(buf(mt,ma,mt,ma,npid,ngid))
    buf=int(h_i,2)
    call proc_output(f,fout,".tmo")
    open(unit=1,file=fout,access="stream",iostat=istat,status="replace")
    if(istat/=0) call echo_red( "*** Error during open: " //trim(fout))
    if(istat/=0) call echo_red("writeHist: error during open: "//trim(f))
    write(1,iostat=istat) buf    
    close(1)    
    if(istat/=0) call echo_red("writeHist: error during write: "//trim(f))
    deallocate(buf)
  end subroutine writeHist

end module oct_pet

